Phase Transitions and Oscillations in a Lattice Prey-Predator Model 
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A coarse grained description of a two-dimensional prey-predator system is given in terms of a 
3-state lattice model containing two control parameters: the spreading rates of preys and predators. 
The properties of the model are investigated by dynamical mean-field approximations and extensive 
numerical simulations. It is shown that the stationary state phase diagram is divided into two 
phases: a pure prey phase and a coexistence phase of preys and predators in which temporal and 
spatial oscillations can be present. The different type of phase transitions occuring at the boundary 
of the prey absorbing phase, as well as the crossover phenomena occuring between the oscillatory 
and non-oscillatory domains of the coexistence phase are studied. The importance of finite size 
effects are discussed and scaling relations between different quantities are established. Finally, 
physical arguments, based on the spatial structure of the model, are given to explain the underlying 
mechanism leading to oscillations. 

PACS numbers: 05.70.Ln, 64.60. Cn, 87.10. +e 



I. INTRODUCTION 

The dynamics of interacting species has attracted a lot 
of attention since the pioneering works of Lotka m and 
Volterra 0. In their independent studies, they showed 
that simple prey-predator models may exhibit limit cy- 
cles during which the populations of both species have 
periodic oscillations in time. However, this behavior de- 
pends strongly on the initial state, and it is not robust to 
the addition of more general non-linearities or to the pres- 
ence of more than two interacting species J3J. In many 
cases the system reaches a simple steady-state. 

A better understanding of the properties of such os- 
cillations is clearly desirable, as such population cycles 
are often observed in ecological systems and the under- 
lying causes remain a long-standing open question pi. 
One of the best documented example concerns the Cana- 
dian lynx population. This population was monitored for 
more than hundred years (starting in 1820) from different 
regions of Canada. It was observed that the population 
oscillates with a period of approximately 10 years and 
that this synchronization was spatially extended over ar- 
eas of several millions of square kilometers ||. Several 
attempts were made to explain these facts (climatic ef- 
fects, relations with the food-web, influence if the solar 
cycle) without success. More recently, Blasius et al. ||] 
introduced a deterministic three level vertical food-chain 
model. The three coupled nonlinear differential equations 
defining the model contain eight free parameters and two 
unknown nonlinear functions. The authors showed that 
an ad-hoc choice of the free parameters and nonlinear 
functions explains the experimental data for the Cana- 
dian lynx. 

In such mean-field type models, it is assumed that the 
populations evolve homogeneously, which is obviously an 
oversimplification. An important question consists in un- 
derstanding the role played by the local environment on 
the dynamics 0. There are many examples in equilib- 
rium and nonequilibrium statistical physics showing that, 



in low enough dimensions, the local aspects (fluctuations) 
play a crucial role and have some dramatic effects on 
the dynamics of the system. Accordingly, a lot of ac- 
tivities have been devoted during the past years to the 
study of extended prey-predator models. The simplest 
spatial generalization are the so called two patches mod- 
els, where the species follow the conventional prey preda- 
tor rules within each patches, and can migrate from one 
patch to the other M. Other works have found that the 
introduction of stochastic dynamics plays an important 
role ||, as well as the use of discreet variables, which 
prevent the population to become vanishingly small. 

These ingredients are included in the so called individ- 
ual based lattice models, for which each lattice site can 
be empty or occupied by one [PHlq] individual of a given 
species or two fl4|,|l5| individuals belonging to different 
species. It was recognized that these models give a bet- 
ter description of the oscillatory behavior than the usual 
Lotka- Volterra (L-V) equations. Indeed, the oscillations 
in these lattice models are stable against small perturba- 
tions of the prey and predator densities, and they do not 
depend on the initial state. It was also found (in two di- 
mensional systems) that the amplitude of the oscillations 
of global quantities decreases with increasing system size, 
while the oscillations persist on local level. It was argued 
that coherent periodic oscillations are absent in large sys- 
tems (although, p do not discard this possibility). In 
fll4| Lipowski et al. state that this is only possible above 
a spatial dimension of 3. In |]ll| Provata et al. emphasize 
that the frequency of the oscillations are stabilized by the 
lattice structure and that it depends on the lattice geom- 
etry. In som e pap ers, the stationary phase diagram was 
also derived |9|,|l5|] , and different phases were observed as 
a function of the model parameters, such as an empty 
phase, a pure prey phase, and an oscillatory region of co- 
existing preys and predators. In |9j, a coexistence region 
without oscillations and a domain of the control parame- 
ter space for which the stationary states depend strongly 
upon the initial condition, were found. 



However, in all the above works no systematic finite 
size studies have been performed, allowing to draw firm 
conclusions on the phase diagram of the models as a func- 
tion of their sizes. It is known p6| , that in ecological 
problems the fact that a system has a finite size is more 
relevant than in most of the cases encountered in statis- 
tical physics, for which one concentrates on the thermo- 
dynamic limit. Particularly, the size dependence of the 
amplitude of the oscillations, as well as a detailed descrip- 
tion of the critical behavior near the phase transitions 
have not been investigated. Another relevant question is 
how much the stationary phase diagrams of these prey- 
predator models have some generic properties or how 
much they depend upon the details of the models. 

The goal of this paper is to study a simple models 
of prey-predators on a two-dimensional lattice for which 
some of the above questions could be answered. Our 
model is based on a coarse-grained description in the 
sense that a given cell models a rather large part of 
a territory and thus can contain many preys or preda- 
tors. Moreover, predators cannot leave without preys 
in a given cell. Those are the main differences between 
our model and Satulovsky and Tome (ST) model ||. 
Nevertheless, it turns out that the stationary state phase 
diagram of the two models are quite different. 

Our model is defined in Sec. O. Although governed 
by only two control parameters, this model exhibits a 
rich phase diagram. Two different phases are observed: 
a pure prey phase, and a coexistence phase of preys and 
predators in which an oscillatory and a non-oscillatory 
region can be distinguished. In some limiting cases the 
model can be mapped onto another well known nonequi- 
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librium model: the contact process (CP) JL7). In Sec 
the properties of our model are analyzed in dynamical 
one and two-points mean-field approximations and no un- 
damped oscillatory behavior is found. In Sec. IV, exten- 



sive Monte-Carlo simulations are performed. It is shown 
that, as a function of the values of the control parame- 
ters, two types of continuous nonequilibrium phase tran- 
sitions towards a prey absorbing state are present. The 
system size dependence of the amplitude of the oscilla- 
tions is studied and several scaling relations between the 
amplitude of the oscillations and the correlation length 
are obtained. In Sec. |V] an underlying mechanism re- 
sponsible for the spatial oscillations is proposed, which 
leads to a qualitative explanation of the properties of the 
phase diagram. In particular, we show that the spatially 
extended aspect of the problem is crucial to have an oscil- 



latory region. Finally, conclusions are drawn in Sec. VI 



II. THE MODEL 



Our system models preys and predators living together 
in a two dimensional territory. This territory is divided 
into square cells, and each of them can contain several 
preys and predators. In this coarse-grained description 



in which each cell represents a rather large territory, one 
can assume that each cell containing some predator will 
also contain some preys. Hence, a three state represen- 
tation is made. Each cell of the two-dimensional square 
lattice (of size L x L, with periodic boundary condition), 
labeled by the index i, can be at time i, in one of the 
three following states: <Xj — 0,1,2. A cell in state 0, 1 
or 2 corresponds respectively to a cell which is empty, 
occupied by preys or simultaneously occupied by preys 
and predators. The dynamics of the system is defined as 
a continuous time Markov process. The transition rates 
for site i are 

i) — v 1 at rate X a (TH,i + n.j ) 2)/4, 

ii) 1 — > 2 at rate A&(nj i 2)/4, 

iii) 2 — ► at rate 1, 

where rii t<7 denotes the number of nearest neighbor sites 
of i which are in the state a. 4 is the coordination number 
of this two dimensional square lattice. 

The first two processes model the spreading of preys 
and predators. The two control parameters, A a and A&, 
characterize a particular prey-predator system. The rea- 
son for considering the sum, n,,i +^2, in the first rule is 
simply that all the neighboring cells of i containing some 
prey, (hence ex, — 1 or 2), will contribute to the prey 
repopulation of cell i. The third process represents the 
local depopulation of a cell due to too greedy predators. 
It can be interpreted as the local extinction of the species 
or as the moving of them to neighboring occupied sites. 
Spontaneous disappearance of a prey state (ctj: 1 — > 0) 
or that of the predators alone (af. 2 — ► 1) is forbidden. 
These assumptions are reasonable because the occurrence 
of these processes is improbable. The rate of the third 
process is chosen to be 1, which sets the time scale. As 
a consequence, A a and Ab are dimensionless quantities. 

The above dynamical rules are an extension of the con- 
tact process model (CP) JL7[ introduced as a descrip- 
tion of epidemic spreading. The CP is a 2-state model, 
<7j = 0, 1; the status and 1 represent respectively the 
healthy and the infected individuals. The CP dynamical 
rules are 

i) — ► 1 at rate A(rij ) i)/4, 

ii) 1 — > at rate 1 . 

An epidemic survives for A > X* CP = 1.6488(1) |L9| and 
disappears for A < X* CP . The transition towards this 
absorbing state is of second order and belongs to the di- 
rected percolation (DP) universality class pq| . 

Our model differs from most of the lattice models pre- 
viously investigated |9H12| by the fact that on each site, 
each species may be represented by several individuals 
rather than just one. In the previously investigated mod- 
els the spreading rate of the preys is simply proportional 
to fii t \. Under this assumption, our model reduces essen- 
tially to the ST model, in which the control parameters 
are defined as c = (1 + A Q + Ab) -1 and p = c(Xb — A a )/2. 



It is worth discussing first the behavior of our model in 
two limiting cases. In the A a — > oo limit the proportion of 
empty cells is negligible since the empty cells are reoccu- 
pied by preys instantly after their extinction. Hence, the 
lattice is completely covered by preys and the a = 2 sites 
behave as the infected species in the CP. Namely, when 
decreasing A& the predator density is decreasing continu- 
ously and vanishes at the CP critical value A^(A a = oo) 
= A£.p. One can think of the A;, — > oo limit in simi- 
lar terms. In this case, the proportion of the prey cells 
(cr = l) should be negligible since the high productivity 
of the predators, while the prey- predator cells should be- 
have as the infected species in the CP. This is indeed the 
case if A a > \* CP , but when A a gets smaller than Xq P , 
the prey density increases again instead of being zero, as 
we shall see later. 



III. MEAN-FIELD ANALYSIS. 

Although apparently simple, there is no way to solve 
analytically the model defined above. However, ana- 
lytic solutions can be obtained by making some approx- 
imations. The simplest one is the one-point mean-field 
approximation in which all spatial fluctuations are ne- 
glected. Thus, the system is characterized by the densi- 
ties of prey, a, and predator, 6, sites 
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which values satisfy the < b < a < 1 conditions by 
definition. In terms of these densities, the mean-field dy- 
namical equations read: 



and 



da 
~dt 



db 



A a a(l — a) — b 



— = \ b b(a — b) - b 
dt 



(2) 



(3) 



Note, that for b = (a = b = 0) initial condition the 
predator (and prey) densities remains 0. 

The (0,0) equations clearly differ from the usual L-V 
ones. The main difference lies in the interaction terms 
as, although a larger prey density increases the predator 
growth rate, the rate of the predated preys only depends 
on the predator density. This is a simple consequence 
of the fact that there are no pure predator sites without 
preys in this model. Thinking of a real prey-predator 
system it makes sense, as a predator has to consume 
a certain amount of preys in a given time to survive, 
independently of the number of preys around it. The 
(1 — a) term in the first equation plays the role of a sim- 
ple Verhlaust factor which assures an upper limit for the 
prey density (a < 1), and similarly the (a — b) term in 
the second equation do not let the density of predators 
exceed that of the preys. 



The stationary states are obtained by setting the left 
hand sides of Eqs. (|| g) to zero. Contrary to the simplest 
L-V equations, qualitatively different stationary states 
are obtained varying the parameters, A a and Ah, as illus- 
trated on Fig. |l|. 



c< 




FIG. 1. Mean-field prediction for the boundary (dashed 
line) between the prey (P) and the coexistence phase (O and 
NO). The dotted lines are the boundaries between the pole 
and node type of stationary state regions. The MC results 
are also depicted for comparison (see Fig. g for the details). 

For < Afe < 1 and A a > 0, the stationary state is a 
pure prey absorbing state a s = 1, b s = 0. For A Q = the 
stationary state is also a prey state, b s = 0, however, the 
value of a s depends upon the initial state. 

In the rest of the plane (A a , Ah), the stationary solution 
is 



(A„ - 1) + y/(A„ - l) 2 + 4A n /Afc 
2A Q 



and 



b s = a s - — 
Ah 



(4) 



(•5) 



which describes a coexistence of preys and predators (co- 
existence phase). 

For Ah ^ 1 the a and b densities are approximately the 
same 



1 



a s = b s + O(-) = { 
Ah 



0(4=) 






for A a > 1 
for A a = 1 
for A a < 1 



(G) 



and as a function of A a they show a mean field CP behav- 
ior as it is expected from the argument given in Sec. |l[ 



In the A a 3> 1 limit (and for Ah > 1) the system is "full 
of preys" , namely 



1 



1 , 1, 



and the predator density reads 
b s = (1 






0(±) 



+ <%) 



(7) 



(8) 



and, as expected, its A b dependence agrees with the pre- 
diction of the mean-field approximation for CP. This 
approximation predicts a second order phase transition 
along the whole A b = 1 line, as in the A b — > 1 limit a and 
b approach linearly the values 1 and respectively 



1 



A„ 



A,, 



1 



Aa + 1 



+ 0((A fc -l) 2 ) 



Ah-1 

A a + 1 



+ 0((A b -l) 2 ) 



(9) 



The behavior of the densities is rather surprising at 
the A a = boundary of the coexistence phase. For 
< X a < 1 and for A b > 1 



i + A "(i + i' " iA " 



(10) 



while the stationary solution, a s , for A a = depends 
on the initial state. Thus the mean field approximation 
predicts a discontinuity of the prey density along this 
boundary. However, the density 6 s = a s — X^ is propor- 
tional to A a and continuous in A a = 0. 

Important quantities are the fluctuations of the prey 
and the predator densities (mean square deviations), 
which are normalized to be size independent for large 
systems 



X P = L 2 ((p - (p)) 2 ), with p = a or b 



(11) 



and () means the time average in the stationary state. 
For A a , A b 3> 1 the majority of the sites are in state 2, 
with a few holes in it, hence one can suppose that the 
holes are independent. Consequently, the number of the 
holes follows a Poisson distribution, from which the av- 
erage hole number equals to the mean square deviation. 
There are L 2 (\ — a) holes made of sites in the state a. L = 
and L 2 (l — b) holes made of sites in the states <7j = or 
1. Thus 



Xa 



1 



and Xb ~ 1 — b s 



(12) 



which is in good agreement with the simulations in a re- 
gion (non-oscillatory part) of the coexistence phase (see 
Fig. I). 

The stability of the stationary state can be analyzed 
by linear stability. One has to investigate the eigenval- 
ues, ei^, of the Jacobian matrix related to the mean field 
equations (J2J, gj) at the stationary densities 



d a a d b a 
d a h d b b 



A„(l- 
A b a s 



2a s 
-1 



-1 
X b a s 



(13) 




FIG. 2. Pole type of approach of the stationary solution 
in mean-field approximation for X a = 1 and A;, = 2, starting 
the system from different initial conditions . Note, that the 
< b < a < 1 conditions are always satisfied. 

It turns out that the real parts of the eigenvalues are 
always negative, assuring the stability of the solutions. 
This mean field approximation do not predict limit cy- 
cles, which would correspond to having an eigenvalue, e, 
with a zero real part. However, in some part of the coex- 
istence phase the imaginary part is nonzero, so the sta- 
tionary solution is approached in spirals (poles), instead 
of straight lines (nodes) (see Fig. |Tj and ||), as it was also 
observed in the ST model ||. Note, that an unexpected 
node region appears for A b > 10. One can consider the 
presence of poles as a hint for the appearance of oscilla- 
tions beyond the mean-field approximation. Notice, that 
in this pole case, the damped oscillations are strong along 
the A b axes (i.e. for A a <C 1). The strength of them can 
be characterized by the ratio of the imaginary and real 
part of the eigenvalues, which has a singularity in the 
A a — v limit. Using (|10|), we obtain 



9(e) 



R(e) 



AX-" 2 ^X 2 -X b + 0{Xl' 2 ) (14) 



for A a <C 1 and A b > 1. In this limit one can derive 
an expression also for the frequency, lj, of the damped 
oscillations 



= |9f(e)| = 2A^ 2 



i-±-+o(xr) 



(15) 



The mean field results can be interpreted in the follow- 
ing way. The approximation predicts two distinct phases: 
the pure prey phase and the coexistence one. It also gives 
some hints for a possible presence of oscillations in some 
parts of the coexistence phase. The phase boundaries of 
the two phases are described by two lines: the A b = 1 



and the A a = 0. Several quantities show a power-law 
behavior close to these boundaries, like b and 1 — a at 
the Xb = 1 boundary, and b, u> and the strength of the 
damped oscillations at the A a = boundary. This im- 
plies that the transitions are of second order, and the 
predator density, b, seems to be a good candidate for the 
order parameter. The order parameter goes to zero at 
the phase boundaries as b ~ (Xb — l) 13 and b ~ Af with a 
mean- field exponent (3 = 1. 

We performed also a pair approximation, in which the 
nearest neighbor correlations are also considered as pa- 
rameters. It turns out that the results differ only quanti- 
tatively from that of the one point approximation. Con- 
trary to U , our system does not show limit cycle behavior 
on the pair approximation level either. 



IV. MONTE-CARLO SIMULATIONS. 



simulations. In one elementary time step one lattice site 
is chosen randomly and its state evolves according to the 
rules defined is Sec. || using rescaled rates (all less than 
1) as transition probabilities. One MC step is defined as 
the time needed such that all the sites have been, on the 
average, visited once. In this paper we always use the 
original time units defined by the model, which can be 
obtained simply by rescaling the time measured in MC 
steps. 

For sufficiently large system, the stationary state does 
not depend on the initial conditions. Usually we filled 
up the lattice completely with preys as an initial state 
and put a few predators on it. To obtain the stationary 
phase diagram and the stationary values of the quantities 
of interest the number of MC steps performed varied in 
the range 10 6 to 10 5 for systems of linear size L = 200 to 
1000 respectively. 



On general grounds, one expects that the fluctuations 
will play an important role in low dimensions. Our model 
is supposed to describe a two dimensional world and 
accordingly, we have performed extensive Monte-Carlo 
(MC) simulations for systems of sizes L x L, L varying 
between 100 and 1000. We used periodic boundary con- 
ditions. 
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FIG. 3. Stationary state phase diagram as obtained by sim- 
ulations. The squares (□) indicates the transition to the prey 
absorbing state (P) for different system sizes (L = 100, 200, 
500 and 1000), and the arrows points to the Xcp value. On 
all figures larger symbols correspond to larger systems. The 
boundary between the oscillatory (O) and the non-oscillatory 
(NO) region of the coexistence phase is determined based on 
Fourier analysis (o) and on the crossover in Xa (•)■ For the 
DP type transition between P and NO, the fitted values of 
Aj(Aa) (x) and the approximation described in Sec. M (solid 
line) is also depicted. 

Although our model is formulated as a continuous-time 
process, an equivalent (at least for not very short times) 
discrete time formulation is more suitable for numerical 
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FIG. 4. Typical behavior of the prey, a, and the predator, 
b (b < a), densities in the oscillatory region of the stationary 
state (X a = 0.8, X b = 100, L = 1000). 

The corresponding phase diagram is depicted on Fig. g. 
Two different phases are present as a function of the two 
control parameters A a and A&: a pure prey phase, a prey 
and predator coexistence phase with an oscillatory and a 
non-oscillatory region. In the oscillatory region, oscilla- 
tions with a well defined frequency were observed in the 
prey and the predator densities (see Fig. ||). Although 
theoretically possible, we never observed an empty lattice 
absorbing state. The reason for that is simply that even 
one surviving prey fill up the system with preys in the 
absence of predators. As Fig. || shows, the locations of 
the different regions of the phase space differ essentially 
from those obtained for the ST model. 

The phase boundaries of the prey phase (see Fig. ||) 
were obtained in the following way. Simulations were 
started at parameter values for which the coexistence is 
maintained practically forever (up to the maximal num- 
ber of MC steps investigated), and we decreased one of 
the parameter values by AA. If the predators were still 



alive after a given time, At, we decreased the parame- 
ter further. The extinction of the predators defines the 
phase boundary. AA was chosen to be in the range 0.005 
to 0.04, while At = 3 x 10 4 MC steps. The result was 
very similar with At — 10 4 and 5 x 10 4 MC steps. The 
definition of the boundary between the oscillatory and 
the non-oscillatory region of the coexistence phase will 
be described later. 



to be the same when the transition line is crossed while 
decreasing A Q at fixed values of A&. The two exponents, 
Pi and 71, are compatible with those obtained for DP in 
2 + 1 dimension [EOj . Thus we conclude, that this absorb- 
ing state phase transition belongs to the DP universality 
class, as expected on general grounds f21|] . Note that DP 
type phase transition in a similar lattice prey-predator 
system has already been observed in 1 dimension |Tq |. 




FIG. 5. The same as on Fig. w but as a function of the 
variables used in the ST model. The triangle represents the 
available part of the phase space. The location of the os- 
cillatory (O) and the non-oscillatory (NO) regions are quite 
different from that of the ST model. 

On Fig. ||, the boundary of the prey phase is displayed 
for different system sizes (L — 100, .., 1000). Apparently, 
in the A a > A& regime the size dependence is negligible, 
but relevant for A a < A&. Note that this strong size de- 
pendence of the boundary coincides with the presence of 
oscillations. 

Decreasing A& at any fixed value of A a , a second order 
phase transition takes place between the coexistence and 
the prey absorbing phases along a transition line A£(A a ). 
As for the mean field case, the predator density is con- 
sidered to be the order parameter. As A& — * AJ(A a ), the 
order parameter, b, and 1 — a go to zero as 



1 - a ~ (A 6 - A£(A a )) ft 



(16) 



As seen on Fig. [|, the values of A£(A a ) obtained by fit- 
ting the data with Eq. [lq are in very good agreement 
with the phase boundary obtained previously. Fitting 
the data leads j3\ « 0.58(1) (with satisfactory precision 
for A a >0.3; see Fig. §). 

In the same limit the fluctuations of the predator 
density also follow a power law behavior, \b ~ (A& — 
A£(A a )) 71 . The exponent has been determined to a good 
precision as 71 ~ 0.35(3) for several values of A a between 
1 and 50. The same behavior has been obtained (only 
for A a = 1 and 3) for the prey fluctuations, Xa, with 
the exponent 71 sw 0.35(5). The critical behavior seems 
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FIG. 6. Prey (open symbols) and the predator (filled sym- 
bols) densities close to the second order phase transition line 
between the prey phase and the non-oscillatory region of the 
coexistence phase. \ a — 0.5(O), 5(0) an d 100(D) while the 
system sizes are L — 200 and 500. The slope of the dashed 
lines is the DP critical exponent j3 R3 0.583. 

For A a — * 0, the transition line, A£(A a ), ends in a spe- 
cial point, (A a = 0, X[ ss 5.0(3)), where all the three 
phases meet. For A > A^, the MC results for the fi- 
nite size dependent phase boundary suggest us that the 
transition happens at A* = in the L — > 00 limit. Ap- 
proaching this transition line, A a — ► 0, the prey density, 
a, does not go to 1 but to a finite value depending on A 
(see Fig^Q). However, according to the results depicted 
on Fig. |8j, the predator density, b, always goes to zero in 
this limit as a power of A a with an exponent fi% s=s 1. Sur- 
prisingly, this second order phase transition to the prey 
absorbing phase does not belong to the DP universal- 
ity class. The presence of power law behavior, however, 
confirms that for an infinite system the transition occurs 
at A* = for A& > A^. This means that for this range 
of A , and for any arbitrary small A a , the coexistence of 
the species is possible providing that the system is large 
enough. 

For A > A^ the fluctuations of the two densities, \ a 
and Xb> behaves similarly. For a given Ab, there is a clear 
crossover at A^(Ab) from a mean field like behavior to a 
regime where the correlations are more important. For 
A Q > A^(Ab) the behavior of Xa and Xb agrees with that 
predicted by mean-field theory, reflecting the fact that in 
this range of A a the dominant behavior comes from the 



noise. As the A a < X°(Xb) condition coincides with the 
presence of oscillations, the crossover point, A^(Ab), is 
taken as the definition of the border between the oscilla- 
tory and the non-oscillatory region. 



not satisfactory enough to obtain the functional form of 
Ki(Xb) (and of the forthcoming Ki(Xb) for i = 2 ,3 and 
4 either). 
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FIG. 7. Prey density, a, as a function of A a for different val- 
ues of A 6 = 3(v), 4 (pentagon), 4.5(0), 5(0), 10(A), 100(D) 
and system sizes L — 200, 500, 1000. The dashed line is the 
density given by the CP. 
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FIG. 8. Predator density, b, as a function of X a for different 
values of X b = 3(x), 5(0), 10(A), 20(O), 50(v), 100(D) and 
system sizes L = 200, 500, 1000 and 2000 only for \ b = 5. 
The A a — > behavior is close to a power law with an exponent 
1 (solid line), while the dashed line is the density given by the 
CP. 

After a proper normalization, the relative fluctuations 
collapse on a single curve for A a < A^(Af,) (see Fig. ||). 
Namely, 



Xb_ 

b 2 



^i(A 6 ) 



Xa 



(1-0)2 



(17) 



where the numerical factor, Ki(Xb), depends only on 
Af,. However, the precision of the simulation results was 



FIG. 9. Normalized fluctuations of the prey (open symbols) 
and predator (filled symbols) densities, Ki(Xb)Xa/(^ — a) 2 and 
Xb/b 2 , which collapse in the oscillatory region. The parame- 
ters are A 6 = 5(0), 10(A), 20(O), 50(v), 100(D) and system 
sizes L = 200, 500, 1000. The dashed lines correspond to the 
mean-field solutions ([_2J). 

The simulation showed that Xp (p = a or &) i s s i zc 
independent as it was expected from its definition ([11]). 
As a consequence, the deviation from the average den- 
sity, a — y/xp~/L pi], is smaller for larger systems and 
evidently scales with 1/L. Certainly, this deviation in- 
creases with the intensity of the oscillations. The above 
finite size behavior is in agreement with the results of 
earlier simulations which claimed that the oscillations 
in the global densities disappear with increasing system 
size |lffl . Our simulations predict more pronounced os- 
cillations for smaller A and for larger Xb- 

The oscillations have to show up also in the correlation 
functions, 

C a {i,T) = ((1 - 5 aj ( t ) t0 )(l - d„ j+i ( t +T),o)) , 

C b (i,T) = {5 ai (t)fi8 a . + .(t +T)a ) , (18) 

where j + i labels a lattice site distant of i lattice spac- 
ing from the site j. C p (p — a or b) depends only on 
i and r due to the homogeneity of the system in space 
and time. For r = the correlation function, C p (i), 
obtained numerically could be fitted by an exponential 
C p (i) ~ exp(— i/£, p ). In the oscillatory region the corre- 
lation lengths of preys and predators differ only through a 
Af, dependent factor, £ a ss ^(Ab)^, and they turned out 
to be proportional to the fluctuations of the prey density, 
£a ~ V^Xa- It means that a more correlated system dis- 
plays stronger oscillations. The reason for that is simply 
that the dynamics of the different sites shows some syn- 
chronization within a correlation length, which results in 
larger oscillations (see Sec. |v| for more details). 



In order to determine the characteristic frequency, 
w p (A a ,Ab), and the amplitude, A p (X a ,Xb) (p — a or 6), 
of the oscillations, we measured the Fourier spectrum of 
the time dependent densities 



S p {uj) 



= lim — 

T^ooT 



T 

E 
t=i 



p(t) exp(iLut) 



(19) 



The presence of oscillations is reflected as a peak at 
nonzero frequency in the Fourier spectrum. Extracting 
this peak from a background noise, enable us to define A p 
and ujp as the zeroth and the first momentum of this dis- 
tribution. This analysis shows clearly that the frequency 
of the oscillations is independent of the system size (see 
Fig. [l(j), and is the same for preys and predators. More- 
over, for a wide range of the parameters in the oscillatory 
phase the frequency, ui = u> a = LJb, is well approximated 
by A a /2. This linear behavior differs from the mean field 
prediction. 
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FIG. 10. Frequency of the oscillations as a function of A a 
and for A b = 4 (pentagon), 5(0), 10(A), 20(O), 100(D) and 
for system sizes L — 200, 500, 1000. For a wide range of the 
parameters the data are close to A a /2 (dashed line). 

In the oscillatory region the oscillations arc present 
for arbitrary large systems, however, their amplitude de- 
creases with increasing system size, as 1/L . At this 
point it is important to emphasize that this fact does not 
imply that only small oscillations are present in large 
systems. Indeed, for a large system the amplitude of the 
oscillations can be made larger by decreasing A Q . On 
the other hand, when increasing A a the amplitude goes 
to zero as a power law which makes difficult to define a 
phase boundary for the oscillations in this way. However, 
there is a simple scaling relation between the amplitude 
and the correlation length for the preys in the oscillatory 
region 

£ « 2 Xa « L 2 A a , (20) 

as it can be observed on Fig. O. The analogous expres- 
sion for the predators is slightly more complicated 



ti 



1 -a 



K 3 (X b ) « X b » K 4 (X a )L 2 A b 



(21) 



with appropriate K 3 (Xb) and K^(X a ) values. 

Another quantity which characterizes the oscillations 
is the time dependent local correlations, C p (t) = C p (i = 
0, r). A similar investigation was made in pdj | with time 
dependent correlations of the average local densities. In 
the oscillatory region C p {t) displays damped, size inde- 
pendent oscillations. More precisely, the time correla- 
tions are size independent for any L > L c (A a ,Ab), while 
for any L < L c the system evolves to the prey absorbing 
state. Clearly, this critical system size is proportional to 
the correlation length, L c ~ £. The size independence 
of C(t) is a simple consequence of the fact that areas 
which arc further than £ apart are uncorrelated. The in- 
vestigation of the time dependent correlations, however, 
provides a rather ambiguous way to define the boundary 
of the oscillatory region. Indeed, one can observe local 
oscillations everywhere in the coexistence phase simply 
because, due to the cyclic dominance nature of the model, 
each site has to evolve in a loop (a = — ► 1 — -> 2 — ► . . .). 
Thus, according to the value of the damping factor, it is 
somehow arbitrary to decide if the state is oscillatory or 
not. 
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FIG. 11. Test of the relation between several characteris- 
tics of the prey population (see Eq. EOh , namely the corre- 
lation length, £ a (A, v)i the fluctuations of the prey density, 
Xa(0, O), and the amplitude of the oscillations, A a (Q, pen- 
tagon) for At, = 5 and 100 respectively. The sizes of the system 
are L = 200, 500 and 1000. 

It is worth noting, that at some particular values of Xb 
(Xb = 10 or 20) and for small A a values (A a < 0.2 or 0.4 
respectively), where the correlation length is comparable 
to the system size (L ~ 500), the system can evolve to a 
stripe like state. In this state 3 stripes of size L, made of 
predator, prey and empty cells, are drifting through the 
system. However, for given A a and Xb values, this behav- 
ior disappears when increasing the size of the system. 



The comparison of the MC results with the mean-field 
prediction shows that the later gives a qualitatively cor- 
rect description of the phase diagram (see Fig. ||), as well 
as of the discontinuity in the prey density, a, along the 
A Q = boundary. 



V. DISCUSSION 

A qualitative understanding of the phase diagram is 
possible. If the birth rates are much larger than the death 
rate (A a 3> 1 and Xb ^> 1) the system is full of preys and 
predators (a « b « 1), while for small values of Xb the 
system evidently reaches the pure prey absorbing state. 

As already discussed in Sec. ||, in the A a — > oo the 
system is full of preys (a — > 1) and the predators behave 
like the infected species in the CP. It means that they 
could survive only for At, > X* CPl where a DP like second 
order transition occurs. This is in agreement with the 
mean field results and with the simulation for A a = 100 
(see Fig. O). 



App + c/A,,. 



parameter the transition occurs at Aj*(A a 
This conjecture is in excellent agreement with the simu- 
lation data for A a > 0.5 with c = 1.28(3) (see Fig. ||) 




FIG. 12. Typical stationary state configuration of preys 
(grey) and predators (black) on a 200 x 200 lattice at A a = 0.9 
and A;, = 100. The white parts represent the empty sites. The 
picture shows the beginning of the invasion of the pure prey 
territory by predators, which were screened by empty sites 
before. 

One can also derive an approximate formula for the po- 
sition of the phase boundary between the non-oscillatory 
and the prey phase A£(A a )- For A a 3> 1, the system is 
almost full of preys (a rs 1) and, in some sense, the dy- 
namics of the predators is close to that of the CP. The 
predators die at rate 1 and spread at rate Xb, but they 
cannot enter into the empty sites. One can introduce an 
effective Xb and describe the process as a CP, namely, the 
predators can enter any neighboring site at this rate. As 
the number of empty sites is proportional to leading order 
to 1/A a , the effective parameter should be Xb = Aj — c/A a , 
where c is a fitting parameter. As this CP displays a 
phase transition at Xb = X^ P , in terms of the original 
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FIG. 13. The same as Fig. |l| but for X a = 0.2 and X b = 5. 
Note, that the predators invade only the fully dense prey areas 
on both figures. 



For A;, ^> 1 the new prey sites are usually immediately 
occupied by predators as well. However, with a small but 
finite probability, a predator site can disappear before the 
predators spread to the new born prey site, and in this 
way, a prey site can be left alone and grow (similarly to 
the Eden model |23]). This rare event is negligible when 
the predator density is large enough and a prey island 
cannot grow for long periods of time. Practically this is 
the case for A a > A^ P . In this case, the number of prey 
sites is negligible small, and the predator sites behave as 
the infected species in the CP. One can see on the Fig. 
and H, that for Xb = 100 the two densities (a rs b) are 
equal to that of the CP if A a > Xq p . However, in the 
vicinity of Xq P the densities are low, which allows an iso- 
lated prey island to grow for a long time. If A a < X^p 
the predator islands are shrinking and, if A a and L are 
not too small, they can survive until a growing prey is- 
land reaches one of them. At this moment, the predators 
invade very quickly the prey territory and increase their 
population size (see Fig. [l2]). These new predator sites 
start to die out leaving a few prey sites alone, and the 
whole procedure starts again. This mechanism insures 
the survival of the predators much bellow the CP critical 
density and results in oscillations in the population sizes. 



,. ; , but not too large, the qualitative picture 



For A fe > A 

is slightly different. As one can observe on Fig. |bj, groups 
of predators are wandering through the system towards 
prey-dense areas. If two fronts of predators meet they 
usually stop moving and the local population of preda- 
tors starts shrinking. The oscillations are maintained in 
a somewhat similar way than for the A& 3> 1 case: these 
predators can only survive if the preys become dense 
around them. This is more probable for larger values of 



A a , and it is also clear that the predators have a better 
chance to survive in larger systems. 

According to the above statements, the key point in the 
underlying mechanism of oscillations is the existence of 
blocked predator islands which are located and trapped 
in sparse prey areas. Indeed, blocked predators in sparse 
prey areas result in growing prey populations; however, 
the resulting dense prey population allows predators to 
move and predate again. This mechanism drives back the 
system to the beginning of the loop. Clearly, predators 
can only be trapped in sparse prey areas if A a is smaller 
or of the order of the death rate, 1. This explains the 
location of the oscillatory region. Note, that the above 
argument is based on the spatial nature of the system, 
suggesting that the spatially extended character is fun- 
damental for the existence of such prey-predator type of 
oscillations. 

This mechanism also provides a qualitative under- 
standing of the key properties of the system. The trapped 
predators can invade the prey area only when the preys 
are dense enough again, which takes a time proportional 
to 1/A a , and leads tow~ A a . According to the simula- 
tions the correlation length, £, increases with decreasing 
A a . Indeed, as A a decreases, the trapped predators have 
to wait longer to escape, hence fewer groups of predators 
survive. This increases the distance between the groups, 
resulting in larger prey islands, which average size is pro- 
portional to £ a - 

When the correlation length is of the order of the sys- 
tem size, there are islands of preys of typical size L, 
extruding the predators out of the system. Hence, the 
condition £ a ~ L characterizes the phase boundary be- 
tween the oscillatory and the prey phase. On the other 
hand, a correlation length of order one (£ a ~ 1), means 
that the noise dominates the system. Thus, £ a ~ 1 char- 
acterizes the boundary between the oscillatory and the 
non-oscillatory region of the coexistence phase. 

As shown by the study of the time dependent corre- 
lations, domains separated by a distance larger than £ a 
oscillate asynchronously around a constant value with the 
same frequency, U)(X a , A&). According to this picture, one 
can derive a more quantitative description for the oscil- 
latory region. Let us assume that for 1< ^ a < L, the 
global densities of each species can be written as the sum 
of local coarse-grained densities at a typical length scale 
£ a . Moreover, we assume that all these local densities os- 
cillate with the same frequency but a different phase, a;. 
In general, the amplitude of the local oscillations should 
depend on the parameters A a and A;,. However, as one 
can observe on Fig. [13 and 13, the predators can only 



a(t) 



2 (f) 2 

2 sin(wi + ai) 
1=1 



(22) 



enter an almost fully dense prey area and the predator 
fronts leave an almost empty field behind them. Hence, 
as suggested by the numerical simulations, everywhere 
in the oscillatory region, the local amplitude for the prey 
density can be considered as a constant, d. Thus 



Supposing that the ai values change much more slowly 
than to, a(t) shows a simple sin behavior for long periods 
of time (see Fig. ||). Thus, for a(t) one can derive the 
value of the density fluctuations, Y a , and the amplitude 
of these oscillations, A a , using ( $A\ ) and (|19|), and take 
the average over all the possible a; configuration taken 
from a flat distribution. This procedure reproduces the 
result of Eq. (GO) up to a multiplicative factor in front of 
the correlation length. 



VI. CONCLUSIONS 

We have studied a two dimensional prey-predator 
model, (size L x L), which exhibits a rich stationary 
state phase diagram. A particular attention has been 
payed to the study of finite size effects, and we were able 
to draw clear cut conclusions concerning the behavior of 
the model both for L finite as well as for the limit L — > oo. 

Three kinds of stationary states can be reached ac- 
cording to the values of the control parameters : a pure 
prey state, and two coexisting prey-predator ones with 
and without oscillations. Two different kinds of second 
order transitions were found when going into the prey 
absorbing phase. The transition between oscillatory and 
non-oscillatory coexistence phase is, in fact, a cross-over 
between two asymptotic regimes characterized by a very 
small and a large correlation length respectively. In the 
oscillatory regime, scaling relations were established be- 
tween several physical quantities. 

A qualitative explanation for the existence of such os- 
cillatory regime is given, pointing out the crucial role 
of the spatial extension of the system. Indeed, the fre- 
quency of the oscillations is determined locally due to the 
dynamics related to blocked predator islands in sparse 
prey areas. Regions of linear size £ a oscillate with the 
same frequency but with different phases. This explains 
the decreasing amplitude of oscillations with increasing 
system size. On the other hand, slowly changing phases 
result periodic oscillations in the overall prey density for 
long periods of time. Moreover, for suitable choices of 
the control parameters one can have synchronized oscil- 
lations with finite amplitude across arbitrary large sys- 
tems. Thus we think that our simple model could offer a 
qualitative explanation for the behavior of the lynx pop- 
ulation problem described in Sec. Q. 
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